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The time to failure, T, of dynamical models of fracture for a hierarchical load-transfer geometry 
CN| , is studied. Using a probabilistic strategy and juxtaposing hierarchical structures of height n, we 

devise an exact method to compute T, for structures of height n + 1. Bounding T, for large n, we 
are able to deduce that the time to failure tends to a non-zero value when n tends to infinity. This 
fj , numerical conclusion is deduced for both power law and exponential breakdown rules. 
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C/3 : I. INTRODUCTION 

Fracture in heterogeneous materials is a complex physical problem for which a definite physical and theoretical 

I ' treatment is still lacking. By "heterogeneous material" we understand a system whose breaking properties (e.g., 

'^ . strengths, lifetimes) depend on time and/or space in a random way 0] . This randomness arises from the many-body 

G ■ interactions among the constituent parts of the system, each one having mechanical properties that can be considered 

independent of -or at least weakly correlated with- the properties of neighboring parts. The term disordered systems 

is also used as a collective name for this kind of material. The presence of disorder alters radically the way the 

rupture process evolves compared to the single-crack growth mechanism operating in homogeneous materials (such 

as glass or alloys). In heterogeneous material (composites, ceramics, rocks, concrete) the process of rupture begins 

with delocalized damage affecting the bulk of the material, and consisting of an enormous number of microcracks 

fT^ ' nucleated at random inside the system. This population of microcracks evolves with time by coalescence and growth 

^D , of individual microcracks until the final rupture point of the system is reached. In the very final stages, the process of 

^O ■ coalescence gives rise to a single (or a few) dominant crack(s) responsible for the macroscopic failure of the material. 

^T^ . The analytical or even complete numerical solution of this complex problem is prohibitive. Nevertheless, our 

Qs ' understanding of fracture in heterogeneous material has improved recently with the development of simple algorithms 

'••>; . to simulate the breaking process. Most of these algorithms are based on percolation theory |^ and include models of 

d ' random resistor networks 0] , spring networks B , and beam networks S . The standard way of solving these models 

H , is through a more or less dense sampling of failure space via Monte Carlo simulation. But fracture in heterogeneous 

i_^ ' materials is, from a statistical viewpoint, a process critically dependent on the tails of the failure distribution and 

these tails are naturally difficult to sample using conventional Monte Carlo methods. It is thus very important to 

develop a set of simple models which can be analyzed either analytically or numerically, with precision and with clear 

^ . asymptotic behavior, in order to guide our understanding of more complex models. 

The load-transfer models belong to this group of simple, stochastic fracture models amenable to either close analyt- 
ical or fast numerical solution, and whose output, spanning many orders of magnitude in sample size, allows a precise 
characterization of the asymptotic behavior. The collective name given to this type of model is fiber-bundle models 
(FBM), because they arose in close connection with the strength of bundles of textile fibers 0, ||§1. Since Daniels' 
and Coleman's seminal works, there has been a long tradition in the use of these simple models to analyze failure of 
heterogeneous materials. 

FBM come in two "fiavors" , static and dynamic. The static versions of FBM simulate the failure of materials 
by quasistatic loading, i.e., by a steady increase in the load over the system up to its macroscopic failure. One of 
the basic outputs is precisely the value of this ultimate strength. Time plays no role in these models, load a is 
the independent variable, and the strength of each element is considered to be an independent identically distributed 
random variable. On the other hand the dynamic FBM simulate failure by stress-rupture, creep-rupture, static-fatigue 
or delayed-rupture, i.e., a (usually) constant load is imposed over the system and the elements break by fatigue after 



a period of time. The time elapsed until the system collapses is the lifetime or time to failure of the set. Time acts as 
the independent variable, and the lifetime of each element is an independent identically distributed random quantity. 

There are three basic ingredients common to all FBM: first, a discrete set of A^ elements located on the sites of a d- 
dimensional lattice; second, a probability distribution for the failure of individual elements; and, third, a load-transfer 
rule which determines how the load carried by a failed element is to be distributed among the surviving elements in 
the set. 

The most common probability distribution function (second ingredient) used to express the breaking properties of 
individual elements is the WeibuU distribution |9|. For the static cases, where the load , a, is the independent variable, 
failure statistics are described by the function P{a) = 1 — exp{ — {a/ao)''}- Here, ao is a reference strength and p is 
the so-called WeibuU index or shape parameter, which in essence controls the variance in the strength thresholds. For 
the dynamic cases, with time as the independent variable, things are more complicated because the failure of each 
element is sensitive to both the elapsed time and its load history. The probability P(i; a{t)) of a single element failing 
at time t after suffering the load history cr(i) is of the form g 

P{t;a{t)) = l-expi-l^ a., [a(r)] drj , (1) 

where Kj{x), j = 1,2, is the hazard rate or breaking rule. To impart to Eq. (p]) the commonly observed WeibuU 
behavior of real materials under constant load a power-law breaking rule is used pO| : 



ki{ct) = 1^0 I— j . (2) 

Here, vq is the hazard rate (number of casualties per unit of time) under the unit load ctq. For constant load, inserting 
Eq. (g) into Eq. (||) gives the WeibuU probability distribution function for the dynamic FBM: 



Pit;a) = l-expl-iyo{ — ]t\. (3) 

The widespread use of WeibuU statistics stems from the experimental fact that real materials follow very closely 
WeibuU probability distribution functions for both the strength and the time to failure of the individual elements m , 

§, 0, iii- 

Besides the power- law breaking rule, Eq. (^, another popular assumption in composites fracture is the exponential 
breaking rule, 



K2{ct) = <t> exp 



a 

V \ — 

CTo 



(4) 



where (p and ry are two positive constants (the amplitude and the characteristic scale of the exponential function). 
This breaking rule has a theoretical support in the apparent necessity of a Boltzmann factor in the load (stress) for 
any thermally activated process (as fracture at the molecular level is interpreted: [0). The substitution of Eq. (Q) 
into Eq. (|l|) does not give a WeibuU function. Nevertheless, the established use of the two breaking rules makes it 
necessary to take both into account, and so we have dedicated Section V to the discussion of the lifetime of sets under 
the exponential breaking rule. 

The most critical of the three basic ingredients of all FBM is the load-transfer rule, where a great deal of the 
physics of the models is hidden. Three end-members are of interest here: the equal load-sharing (ELS) rule, the local 
load-sharing (LLS) rule, and the hierarchical load-sharing (HLS) rule. In the ELS rule, which can be thought of as 
a mean-field approximation, the load supported by failing elements is shared equally among all surviving elements. 
In the LLS rule, the load of failing elements is accommodated by a neighborhood whose exact definition depends on 
the geometry and dimensionality of the underlying lattice. In the HLS rule the scheme of load transfers follows the 
branches of a fractal (Cayley) tree with a constant coordination number. Common to all three load-transfer modalities 
is the fact that broken elements carry no load. 

ELS models have been used to predict failure under tension in elastic yarns and cables with little or no twist, because 
in these arrangements the load supported by a failing fiber or cable strand is shared equally by all the remaining 
fibers or strands in the bundle. The conjunction of a loose arrangement and a load under tension facilitates this 
global-range load redistribution scheme. 

LLS models have their natural field of application in the failure of composite materials, and more specifically in 
fiber-reinforced composites with brittle fibers embedded in a stiff matrix p3|. There, as fiber breaks appear, the 



matrix serves the important function of transferring the shear traction generated in the matrix at the point of a fiber 
break to the neighboring fibers, with most of the load going to the nearest neighbors. This arrangement results in a 
very short-range load redistribution, both laterally across fibers and longitudinally along the fiber axis. 

More important from a geophysical point of view and for this paper is the HLS rule recently introduced by Turcotte 
and collaborators in the seismological literature llj] . In this load-transfer modality the scale invariance of the fracture 
process is directly taken into account by means of a hierarchical load-transfer scheme following the branches of a 
fractal tree. An important property of the HLS scheme is that the zone of stress transfer is equal in size to the zone 
of failure, and this nicely simulates the Green's function associated with the elastic distribution of stress adjacent 
to a rupture |^. The fractal tree structure used to redistribute loads is a mere construction useful to envisage 
the way loads from breaking elements are transferred to unbroken elements. A basic aspect of the topology of the 
hierarchical structure is the number of elements directly linked together; this defines the coordination, c, of the tree. 
That is, c fibers could be assembled to form a bundle which would behave as if it were itself a fiber. Then, c of these 
second-generation fibers could themselves be assembled to form a bundle which would act as a third-generation fiber, 
and so on. This hierarchical assemblage can be continued indefinitely and an index n is used to describe the level 
within the hierarchy or, equivalently, the height in the tree structure. So n = refers to the individual elements in 
the system, n = 1 refers to the first level in the tree, etc. For c = 2, n = implies individual elements, n — 1 implies 
pairs of elements, n = 2 pairs of pairs of elements, etc. Thus, an nth order tree with coordination number c would 
contain iV = c" elements. 

Although we would stress below the importance of the analytical solution of the FBM, it is enlightening to show 
how these models can be solved using Monte Carlo techniques because that would facilitate the understanding of 
parts to follow. We will focus on the dynamic FBM, as this is the type of problem we want to address in this paper. 
Consider a set of N elements arranged on the sites of a lattice. The general Monte Carlo recipe goes as follows: (1) 
Assign random lifetimes ti to the i = 1, ..., N individual elements, as drawn from Eq. (ffl) under unit load; (2) Advance 
time an amount equal to the lifetime of the shortest- lived non-failed element in the set, say 6; (3) Reduce lifetimes 
of all remaining elements by an amount K{ai)S, where k{x) is either the power-law or the exponential breaking rule; 
(4) Transfer load from the failing element to other sound elements in the set according to a preset load-transfer rule 
(ELS, LLS, or HLS); (5) Proceed to step 2 if at least one element is unbroken, or end if the system has collapsed; (6) 
Add together all the individual Ss to obtain the time to failure T of that particular realization of the system. This 
way of applying the Monte Carlo method is what we will refer to as the standard method. 

Among the different results that one can obtain from the analytical or numerical solution of the fiber-bundle models 
(for a review, see B] ), here we are mainly interested in the asymptotic strength (static FBM) and asymptotic time 
to failure (dynamic FBM) of the system. The asymptotic strength is defined as the maximum load that an infinite 
system can support before all its elements break. The asymptotic time to failure or lifetime is the minimum time one 
has to wait until an infinite system collapses by all its elements breaking by fatigue. These are themselves important 
questions from an engineering point of view. Table 1 gathers the main asymptotic results for the different FBM, 
including the dynamic HLS model, to which this paper is dedicated. It has been known since the work of Daniels 
|[7| that the static ELS fiber-bundle model has a critical point in the sense that for an infinite system there is a zero 
probability of breaking the system when applying a load a less than a critical value CTc, and a probability equal to one 
to break the system if the applied load is bigger than the critical load. This is valid for any probability distribution 
function satisfying some very mild conditions [|7|. The critical strength CTc quoted in Table 1 is only valid, however, for 
a WeibuU function. As for the dynamic ELS model, Coleman |8| proved rigorously a comparable result, namely that 
there exists a critical time Tc below which an infinite system under dynamic rules has a zero probability of collapsing 
and above which the system collapses with a probability of one. Tc varies with the assumed probability distribution 
function, but otherwise the critical-point result is independent of it (the value quoted in Table 1 is for a power-law 
breaking rule) . Smith and Phoenix |16| give a summary of the main asymptotic results for the static and dynamic 
ELS FBM. Regarding the asymptotic properties of the LLS models, the work of Smith and co-workers, and of Harlow, 
Phoenix and co-workers has been fundamental. They proved that neither the static nor the dynamic LLS FBM have 
a critical point. For the static case, the strength of the system goes to zero as the size of the system is increased. More 
specifically, Cc oc 1/lnN (conjectured by Harlow, Phoenix and Smith |I^; proved by Harlow llq]). For the dynamic 
case a similar result holds (conjectured by Tierney |19]; proved by Kuo and Phoenix CQ]). Sec |2^| for a review of 
the static LLS models. 

The static HLS model was shown to lack a critical point by Newman and Gabrielov |23] . In this case the reduction to 
zero of strength in relation to system size is very slow, CTc oc 1/lnlnN, but strictly speaking the strength of an infinite 
system is zero. The dynamic HLS model was introduced in the geophysical literature in reference p3| . Afterwards, 
Newman et al p5| used this dynamic HLS model with the specific aim of finding out if the chain of partial failure 
events preceding the total failure of the set resemble a log-periodic sequence. This was motivated by the amazing fit 



of this type obtained in Ref. |^^ to data of the cumulative Benioff strain released in magnitude > 5 earthquakes in 
the San Francisco Bay area before the October 17, 1989 Loma Prieta earthquake. In the analysis of ||l^, it appeared 
that, contrary to the static model, the dynamic HLS model seemed to have a non-zero time to failure as the size of 
the system tends to infinity. This issue was also studied in [ p5[ by using a renormalization approach. This behavior 
seems odd because, as Table 1 shows, there is a symmetry, for a specific load-transfer rule, between the static and 
the dynamic cases: the ELS models have a critical point both for the static and the dynamic cases; the LLS models 
have no critical point behavior either for the static or the dynamic case. The static HLS model has no critical point, 
and so it would seem natural that the dynamic HLS would not have a critical point either. Here we present an exact 
iterative method to compute the time to failure of sets of elements with a hierarchical modality of load-transfer. (A 
preliminary account of the behavior of the dynamic HLS model under the power-law breaking rule has been recently 
published p6[). Due to the fact that the exact method is too time consuming to yield useful asymptotic results, we 
also present here rigorous upper and lower bounds to the lifetime of large dynamic HLS sets. From the behavior of 
the lower bound we conclude that the dynamic HLS model has indeed a critical point, that is, its time to failure is 
non-zero for an infinite system. 

This paper is organized as follows. In Section II we review the continuous formulation of the ELS model given 
by Coleman g]. This is useful for understanding the probabilistic approach used in the rest of the paper. This 
approach is compared with the standard method in a Monte Carlo simulation, to illustrate that both are equivalent. 
Section HI contains the core of the iterative method to exactly calculate the time to failure of dynamical HLS sets. 
The strategy of the juxtaposition of configurations is explained and the need of defining "replica" configurations is 
introduced. This leads to the concept of primary diagram from which the value of the 6s and of T, for a given n, 
are exactly calculated. From a primary diagram one obtains an easier one called reduced diagram, which is used to 
build the primary diagram of the next level n -\- 1. In Section IV, aiming at simplification, we introduce the concept 
of effective diagram, as the averaged form of a primary. In these diagrams each stage of breaking is represented by a 
unique effective configuration whose decay width is obtained by an appropriate average of the various decay widths 
existing in the primary. Depending on the type of means used, upper or lower bounds for the T of the next height 
are obtained. So far the power law breaking rule is used in the quantitative calculations. Section V is devoted to 
the exponential breakdown rule specifics. Two Appendices have been added: In Appendix A we detail the number 
of replicas demanded for a general coordination, c; in Appendix B, we show the reason why the three types of means 
used in the averaging effectively work to provide rigorous bounds. 

II. THE ELS MODEL. THE PROBABILISTIC APPROACH 

In Ref. B , in the context of his statistical theory for the time-dependence of mechanical breakdown in bundles of 
fibers at constant total load, Coleman defines an "ideal bundle" as one fulfilling precisely the same premises that the 
ELS model presented in Section I. Following the work of this author, let us call N the size of the bundle (or set) 
at i = and n(t) the number of filaments (or elements) which survive to t without breaking; the lifetime, T, of the 
bundle is defined as the time required for n{t) to reach zero. We will call CTq the fixed load attributed to every single 
element at t — 0. The hypothesis of the ELS model implies that the actual load in a particular unbroken filament at 
time t is 

(7 = — T^. (5) 

n{t) ^ ^ 

Thus the lifetime of a very large ideal bundle formed by fibers of equal length, Z, may be calculated from 

dn , / s ,„\ 

- -^ ="''«('^): (6) 

where ^(cr) is a phenomenological function. The product I k is called the hazard rate. In polymeric fibers, k, can be 
satisfactorily represented by an exponential function of a : 

K{a) = — , (7) 

a 

where a and /3 are parameters that determine the behavior of the fibers under any load a. Eq. (R) represents the so 
called exponential breakdown rule. An alternative also widely used is the so called power-law breakdown rule: 



Ha) = -(-) . (8) 

a \(7oJ 

From (P), (p) and (M) with x — {j3 Gq N)/n, and the conditions n(0) = iV, n(T) = 0, one deduces 

„ a f°° dx e^^ a ^ , „ X 

r=T / ^-jE,i-(3ao). (9) 



Using (H) instead of (^ one obtains 



(a/0 . X 

T = -5-^-^ . (10) 



Eq. (P) is similar to a radioactivity equation in which I ■ k stands for the decay rate, F, of one nucleus. In the ELS 
case it is not of much interest to lose this elegant continuous formulation. However for other load-transfer schemes, 
such as the HLS, this analogy with radioactivity is useful, but similar continuous differential equations can not be 
formulated anymore. Thus, the discrete version of this probabilistic philosophy applicable to any load transfer scheme 
was developed in Ref. |27] and will be used throughout this paper. It represents an alternative to what we have called 
the standard method |15| commented in Section I, in which the random thresholds for breaking are assigned at the 
beginning and the process of breaking is henceforth completely deterministic. Both points of view are equivalent. In 
the following we will use non-dimensional magnitudes. In particular 

(a/l) = 1, Co = 1, and /3 ao = I ■ (11) 

Note that this would be equivalent to adopting, with the notation of Section 1, h'o = ag = (j) = t] = 1. 
In the probabilistic approach p7|], in each time increment, defined as 

one element of the sample decays. The index j runs along all the surviving elements. Using Eqs. (|^), (§) and ( [ll|) we 
have 

Tj=<jP {orTj=e''^) . (13) 

The probability of the specific element, m, to fail is 

Pm=T^-S . (14) 

Equation (n3) is the ordinary link between the mean time interval for one element to decay in a radioactive sample 
and the total decay width of the sample. The time to failure, T, of a bundle (set of elements) is the sum of the N 6s. 
It is instructive to apply Eq. (n3) to the ELS case where Tj does not depend on j because every surviving element 
bears the same load. Here, in the fcth time step, the number of survivors is Uk ~ N — k and the individual load is 
Cfe = N/{N — k). Then for the power law rule 

1 {N - k)P (N - k)P-^ 

Ok = - 



with fc = 0, ..., iV- 1 and 



N -k NP NP 



N-l 



T=J26, (15) 



fc=0 



li p — 2, T — i(l + jr). For a general value of p, using Stolz's theorem [ESl we find 



lim T = ^ , ' ^ > - , (16) 

iv^oo NP-{N-l)P p' ^ ' 



which coincides with Eq. (0 



When dealing with the exponential breakdown rule, proceeding analogously one easily checks that the sum of the 
series of Ss, for sufficiently high N, provides the same result as the exponential integral of Eq. (^). 

It is also instructive to compare the results obtained from Monte Carlo simulations in the calculus of T in two 
ways: (a) by using the standard procedure, i.e., of assigning random individual lifetimes at the beginning of each 
simulation and proceeding deterministically; or (b) by using a probabilistic point of view, i.e., from Eq. (n2h and 
Eq. (UJ). This comparison is shown in Fig. H^ for HLS sets of A^ = 128 and N = 512 elements (with c — 2, p ^ 2). 
Note the significant reduction in the dispersion thickness obtained by using this second method. This contrast tends 
to decrease for growing N and growing p. Obviously, the Monte Carlo strategy can be applied for any modality of load 
transfer in the framework of the probabilistic method. The inconvenience lies in the very essence of these simulations, 
i.e., their moderate accuracy and large cost for large sets. In this paper, we will show how to apply the probabilistic 
method to the HLS transfer modality, in order to obtain an exact algebraic method for the lifetimes, and how to 
explore the asymptotic values of T when N tends to infinity. 

III. AN EXACT ITERATIVE METHOD FOR HLS DYNAMICAL MODELS 

To give a perspective of what is going on in the rupture process of a hierarchical set we have drawn in Fig. g 
the three smallest cases for trees of coordination c — 2. Denoting by n the number of levels, or height of the tree, 
i.e. N = 2", we have considered n = 0, 1 and 2. The integers within parenthesis (r) account for the number of 
failures existing in the tree. When there are several non-equivalent configurations corresponding to a given r, they 
are labeled as (r, s), i.e, we add a new index s. The total load is conserved except at the end, when the tree collapses. 
Referring to the high symmetry of loaded fractal trees, note that each of the configurations explicitly drawn in Fig. S 
represents all those that can be brought to coincidence by the permutation of two legs joined at an apex, at any level 
in the height hierarchy. Hence we call them non-equivalent configurations or merely configurations. In general, each 
configuration (r,s) is characterized by its probability p(r^s), '^gPirj-s) = 1, and its decay width r(r, s). The time 
step for one-element breaking at the stage r, is given by 

^r-j:Pir,s)^^. (17) 

This is the necessary generalization of Eq. (n2[) due to the appearance, for the same r, of non-equivalent configurations 
during the decay process of the tree. In cases of branching, the probability that a configuration chooses a specific 
direction is equal to the ratio between the partial decay width in that direction and the total width of the parent 
configuration. And the probability of a given configuration p(r, s) is given by the sum, extended to all its possible 
parents, of the product of the probability of each parent times the probability of choosing that specific direction. 

We will compute at a glance the Ss of Fig. in order to analyze the general case later. To be specific, we will always 
use p = 2. For n = 0, we have r(0) ^ l^ and Sq = ^ = 1 ^ T. For n = 1, r(0) = l^ + l^ = 2, So ^ h ^(1) = 2^ 
Si = i; and hence T = i + i = |. For n = 2, r(0) = 1^ + 1^ + 1^ + 1^ = A, So = j; T{1) = 2^ + l^ + l^ ^ 6, Si ^ i. 
Now we face a branching; the probability of the transition (1) -^ (2, 1) is | and the probability of the transition 
(1) -^ (2, 2) is §; on the other hand r(2, 1) = 2^ -h 2^ = 8 = r(2, 2), hence ,52 = |-| + |-| = |- Finally ^3 = 1^ and 
the addition of Ss gives T = ||. 

Now we define the replica of a configuration belonging to a given n, as the same configuration but with the loads 
doubled (this is because we are using c = 2). The replica of a given configuration will be recognized by a prime sign. 
In other words (r, s)' is the replica of (r, s). Note that when a configuration represents the state of complete collapse, 
it and its replica are the same thing. When dealing with the power law breakdown rule, any decay width, partial or 
total, related to (r, s)' is automatically obtained by multiplying the corresponding value of (r, s) by the common factor 
cP = 2P = 4. This also implies that p{r, s) = p{r, s)' . In the exponential rule, this does not work and the widths of the 
replicas have to be specifically calculated (this is explained in Section V). The need to define the replicas stems from 
the observation that any configuration appearing in a stage of breaking r of a given n, is built as the juxtaposition 
of two configurations of the level n — 1, including also the replicas of the level n — 1 as ingredients of the game. In 
Fig. g, one can observe the explicit structure of the configurations of A^ = 4 (or of iV = 2) as a juxtaposition of those 
of TV = 2 (or of iV = 1) and its replicas. From this perspective, we notice that the total number of configurations 
appearing in the fracture process of a tree of height n (omitting the totally collapsed one), J\fn, is equal to 



In this formula the first term represents all the possible combinations (with repetition) of pairs of ordinary configu- 
rations of the height n — 1. The second term represents the configurations formed by juxtaposing a collapsed tree of 
height n — 1 together with any of the Mn-i replicas of the previous height. Thus, 

M. . •^"-^(^"-^ + ^^ . (18) 

Feeding TVq = 1 into Eq. (|l|), we obtain M = 2, M = 5, A/'a = 20, M = 230, TVs = 26795, We ~ 3.59 x 10^ 
TVt ~ 6.45 X 10^^, etc. It is clear that the amount of configurations to deal with soon constitutes an unsurmountable 
problem. 

The single-element breaking transitions in configurations of height n can be only of three types. Type a transitions 
correspond to the breaking of one element in a half of the tree while the other half remains as an unaffected spectator. 
Type b transitions correspond to the decay of the last surviving element in one half of the tree, which provokes its 
collapse and the corresponding doubling of the load borne by the other half. In these two cases, the transition width 
coincides with that already obtained when solving the level n — 1. Finally, type c transitions correspond to the scenario 
in which one half of the tree has already collapsed and in the other half one breaking occurs. In this third case the 
decay width is that of a replica of the level n — 1, which as said before, is a common factor c'' times the ordinary 
width. This holds for any height n and allows the computation of all the partial decay widths in a tree of height n 
from those obtained in the height n—1. This is illustrated in Fig. 0, for the three types of transitions for trees with 
n = 3 {p — 2 has been used). 

In Fig. and Fig. we have drawn the different configurations of low n explicitly, that is, by representing them 
as small fractal trees at different stages of damage. It is convenient, for reasons of economy, to introduce a symbolic 
notation for the configurations so that the complete process of breaking of a tree of height n adopts a more compact 
look. This is shown in Fig. 0. There the different configurations of n = 3, are labeled by the integers within the boxes. 
The two parentheses at their right, with their respective integers, represent the two n — 2 juxtaposed configurations 
forming that of level n = 3. This information of the previous height will be called the genealogy. Time is assumed 
to flow downwards. The numbers accompanying an arrow connecting two boxes stand for the decay width of that 
transition. Coming bac k to Fig. H one recog nizes ther e th at tho se e xplici t transitions are nothing else but what in 



Fig. ^ is represented as 3,1 — > 4,2 , 3,1 — > 4,1 and 4,1 — > 5,1 . A diagram like that of Fig. |j is called a 



primary, because it is formed by the juxtaposition of all possible configurations of the previous height. Thus, each 
configuration belonging to a primary diagram has a specified genealogy. This allows the computation of all the decay 
widths of the diagram. As foreseen in Eq. (p^), not counting the totally collapsed configuration, Afs = 20. The sum 
of all the partial widths of a parent configuration in a branching is always equal to the total decay width, F, of the 
parent. From this primary width diagram one deduces the probability of any primary configuration at any stage r of 
breaking, and consequently Sr is obtained using Eq. ([l7|). Finally, by adding all the 5s we calculate T{n — 3). 

After a primary diagram has been obtained, i.e., after calculating all its decay widths, it can be simplified. The 
idea is to fuse, at each r, all the configurations having the same total decay rate, F. Once fused, these configurations 
have a probability equal to the sum of the old probabilities, and obviously maintain the same F. A primary diagram 
simplified in this way will be called a reduced diagram. An element of a reduced diagram resulting from a fusion has 
no genealogy in the sense that it does not derive from one but from several juxtapositions. The genealogy was used 
in the calculation of the primary diagram. The later fusion does not require any other independent information. To 
illustrate the concept of what a reduced diagram is, let us look again at Fig. g. For n — and n = 1, for each r there 
is only one configuration and hence primary and reduced diagrams are identical. For n — 2, for r = 2 there are two 
configurations in the primary diagram, but they have the same width, specifically, for p — 2, r(2, 1) = F(2, 2) = 8. 
Thus these two configurations can be fused and the resulting effective diagram is a chain of five elements, i.e., the 
branching disappears. Performing this task with the n = 3 of Fig. H one would obtain the reduced diagram of Fig. ra. 
The total number of configurations appearing in the reduced diagrams, A/"^, does not derive from a closed formula as 
occurs with A/'„. However, it can easily be derived by means of the computer and obtain: Af{ = 1, M2 = 2, Af^ = 10, 
Afi = 36, A/'s = 202, A/'g = 1669, Mj = 16408. The important point is that one can use a reduced diagram of the 
level n to build a primary of the level n + 1 obtaining the exact information of the new level. After calculating that 
primary, by fusing again configurations of equal F, one would obtain the reduced diagram of the height n + 1. 

By iterating this procedure, that is by forming the primary diagram of the n + 1 height by juxtaposing the 
configurations of the reduced diagram of the height n, we can, in principle, exactly obtain the total time to failure 
of trees of successively doubled size. In spite of the great simplification obtained when using reduced diagrams, the 
problem of dealing with a vast amount of configurations still remains. This fact eventually blocks the possibility of 
obtaining exact results for trees high enough as to be able to gauge the asymptotic behavior of T in HLS sets. A 



few examples of exact results, for c = 2, p = 2 are T{n = 3) = ^^, T{n = 4) = 463oo977698976000 "^(^ = 5) = 
0.420823219104814. The iterative procedure was programed in Mathematica 3.0 with infinite precision and took 10 
minutes CPU time for n = 5. 

IV. BOUNDS FOR THE TIME TO FAILURE OF THE HLS MODELS 

As seen above, whenever in a primary diagram one fuses configurations of the same F, no information is lost and the 
calculation of the time to failure remains exact. In spite of this simplification the magnitude of the Bayesian problem 
becomes huge even when dealing with a moderate n. That is why we have looked for alternative approximation 
procedures to estimate T. In fact the most important goal, as explained in Section I, is to find out if the T of very 
large HLS sets tends to zero or, on the contrary, remains finite. With these points in mind, we have found that a 
drastic but appropriate simplification of the primary diagrams, in which one averages all the configurations of a given 
r into a unique configuration with an effective decay width, leads to obtaining, in the subsequent heights, values of 
T systematically lower (or higher), than the exact result. As this fusion leads to only one effective configuration, it 
will have probability 1. The value of its decay width will be called a^. Such "chain" diagrams will be called effective 
diagrams. For n = 3, this is drawn in Fig. 0. These effective diagrams, which substitute the previously defined 
reduced diagrams, are used exactly in the same way, i.e., to calculate a new (approximate) primary diagram of the 
next height. The economy obtained by using effective diagrams is obvious. As the number of effective configurations 
of a level n — 1 is 2"~^, the primary of height n, built from this effective diagram of height n — 1, will have a number 
of configurations TV", given by 

o2n — 2 _i o on— 1 

M" = ±f^^— ; (19) 

that is A^i" = 2, AA^' = 5, 7V^' = 14, TV^' = 44, Mi; ^ 152, Af^' ^ 560, TVf = 2144, etc. 

It is clear that for n = 0, 1 and 2, the reduced diagrams and the effective diagrams are identical, i.e., a^ = F(r). 
The point is to define a^ for n > 3 so that the T(n > 4) are lower (or higher) than its exact result. 

A trivial option is to define 

ar ^Tmaxir) Or (Tminir)) (20) 

i.e. by assuming that the only configurations formed during the breaking of the tree are those of the maximum 
(minimum) value of F. As it is easy to foresee, the use of Eq. (|2^) leads to poor bounds. In fact the lower bound 
goes quickly to zero. We have found that good lower bounds are obtained by using effective diagrams where a^ is the 
arithmetic mean (AM), 

ar{AM) = Y,Pir,s)T{r,s) , (21) 

S 

or even better by using the geometric mean (GM)^ 

a,(GM) = nr(r,s)P('^^^) . (22) 

S 

Good higher bounds are obtained from the harmonic mean (HM)^ 

aAHM) . y '^ . (23) 

In Appendix B, we analyze why bounds result. Note that given the primary diagram of a height n, which leads 
to c" (5s, the elements forming the effective diagram defined with the aim of obtaining higher bounds are exactly 
Ur — l/Sr- The fact that the ar{i), i = AM,GM, and HM lead to bounds in the form explained above is in 
qualitative concordance with the inequality 

T„iin{r) < ar{HM) < ar{GM) < ariAM) < T,naxir) , (24) 

which is always a mathematical fact. The bounds obtained from these formulae for c = 2, p = 2 are plotted in Fig. M, 
together with points representing Monte Carlo results. As the bounds based on the GM and on the HM are the most 



stringent, they will be called Ti and Th respectively. The detailed behavior of Ti has been analyzed in a log- normal 
plot of Ti — Ti^oo against the number n of levels of the tree. Tj^oo is a constant obtained from a fit of the data points to 
the exponential function ae"'''^"^"''-' shifted downwards by this amount Ti^oo (o, i> and no are three fitting parameters 
of no interest here) . We have performed a careful sensitivity analysis of the four-parameter exponential fitting because 
the success of this exponential decay to a non-zero limit is the hallmark of the claim. Table B records Ti^aa obtained 
from an exponential fit to the last k data points. The first T; oo column is for a fit using up to a maximum level of 
n — 2Q (hence the notation Umax in the table); the second Ti^oo column is for the same fit but dropping the n — 20 
value; for the third TJ co column we have also dropped the n = 19 data point. It is clear from the trend in the three 
T; oo columns that a saturation towards Ti^o — 0.32537 ± 0.00001 occurs when using only information of big trees 
to perform the non-linear fitting. A similar analysis of T/j, leads to Th^oo — 0.33984 ± 0.00001. The quality of this 
exponential fit is also shown in Fig. pi Similar fittings of the Monte Carlo data points are inconclusive, due to the 
intrinsic noisiness of the MC results and the limited size of the simulated sets [N < 2^^ elements). What this result 
implies is that a system with a hierarchical scheme of load transfer and a power-law breaking rule (c = 2, p = 2) 
has a time to failure for sets of infinite size. Too, such that 0.32537 < Too < 0.33984. Thus, there is an associated 
zero-probability of failing for T < Too and a probability equal to one of failing for T > Too • The critical point behavior 
is thus numerically confirmed. 

V. EXPONENTIAL BREAKING RULE 

When dealing with the exponential breaking rule in the probabilistic approach, one has to use Eq. (R) for the 
hazard rate function. For the specific value of the parameters as fixed in Eq. (pT|), we have 

Vj = e"^ . (25) 

The problem with this breaking rule is that the values of the decay widths appearing in a diagram where the loads are 
doubled, i.e., in a diagram replica, are not obtained by multiplying the normal ones by a fixed constant, as occurred 
with the power-law breaking rule. This is easily checked in Fig. where the primary diagram for n — 2 and its replica 
are shown. The notation is equal to that of Section II. The values of the decay widths here are dictated by Eq. (Eq). 
Several comments are in order. We see that the structure of the diagrams is equal to those appearing for the power 
law case because this is independent of the breaking rule assumed. In Fig. M (b) we have drawn explicitly a replica, 
that is a diagram in which the system instead of starting with individual loads CTq = 1 , starts with doubled individual 
loads a'^ = 2ao — 2. We also see that, just in the same way as the quantitative calculation of the primary of Fig. pi (a) 
required the knowledge of the information of the previous height and of its replicas, the calculation of the diagram of 
Fig. |9| (b) requires the knowledge of the primed elements and the double-primed elements (i.e. with loads multiplied 
by 4) of the previous height. Thus, suppose that we want to calculate T up to the height n = A. This demands the 
knowledge of the reduced diagram of n = 3 and of its replica. We will denote them by {3} and {3}'. To obtain {3} 
we need to know {2} and {2}', and to obtain {3}' we need to know {2}' and {2}". Going backwards up to n = 0, we 
observe that the scheme of information needed looks like the following triangular array : 



{0} = 


= e {0}' = 


= e' {0}". 


.e^ {0}"' = e« 


{1} 


{1}' 


{1}" 


{1}'" 


{2} 


{2}' 


{2}" 




{3} 


{3}' 






{4} 









{0}"" = e 
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In other words, to obtain the T of a given height n we have to explicitly calculate the primary diagrams of the 
previous values of n, starting from n = 0, up to a loading n times the usual diagram with CTq = 1. The primary 
diagrams of low n are calculated at once, hence the extra work with respect to the power law case is not too much. 
When explaining the above mentioned triangular array we were referring to exact primary diagrams, taking for granted 
that our aim was to obtain exact results. Obviously this changes if our aim is the calculation of bounds; then one 
would proceed by averaging widths for each r, that is by calculating means and dealing with effective diagrams. 

In Fig. nol we have drawn the lifetimes, T, for trees of height n. The circles are results obtained from Monte Carlo 
simulations. For the exponential breaking rule, the lower bound based on the arithmetic mean, curve (3) goes to zero. 
Thus the only lower bound that remains useful is that based on the geometric mean. Again, it will be called T;. By 
fitting the data T; by an exponential function of the form T; = Tj^oo + ae~^'^"~"°'' we observe a clean saturation of 
the asymptotic time to failure towards Tj.oo = 0.05285 =F 0.00001; analogously, we obtain Th,oo = 0.08825 T 0.00001. 
Hence, the critical point behaviour is also numerically confirmed for the exponential breaking rule. 



VI. CONCLUSIONS 

In this paper the time to failure, T, of hierarchical load-transfer models of fracture has been studied. Initially 
we have explained in detail the so called probabilistic approach to load-transfer dynamical models as opposed to the 
standard approach, in which random lifetimes are assigned to the elements of the set and the process of fracture evolves 
dcterministically. We have emphasized that when viewed from the probabilistic point of view, the calculation of T is 
analogous to the computation of the total decay time of a radioactive sample. In fact the terminology of radioactivity 
appears throughout this paper. We have shown that the calculation of T, using Monte Carlo simulations, has a 
smaller dispersion if one adopts the probabilistic approach. 

Then, we have devised an exact method to compute T of hierarchical structures of size N — c". The number of 
elements of the set is N, c is the coordination of the tree and n the height of the fractal tree. The method is iterative, 
i.e., for a given c, it allows the computation of a tree of height ri + 1 once one has calculated a tree of height n. In this 
context, the sentence "a tree is calculated" means that one knows the value of all the partial decay widths between all 
possible configurations appearing during the breaking process of that tree. Once this information is known, one easily 
calculates the probability of reaching each configuration, and the individual values of each S, i.e., the one-element 
breaking time. The key of the method derives from the observation that the structure of the configurations of the 
n + 1 type are a mere juxtaposition of c configurations of the n type. In this juxtaposition, the so called replicas 
also play a role. The quantitative information of how replica configurations behave is explained for the two relevant 
breaking rules: the power law and the exponential. In the power law breaking rule any decay width of a replica is 
just a common factor times the original value. In the exponential breaking rule, on the contrary, the decay widths of 
replicas have to be individually calculated. 

The iterative process, including the information of the replicas, can be easily processed by a computer. It allows 
the exact calculation of T for moderate heights n. An exact simplification, denoted as reduction, is introduced to 
diminish the magnitude of the information to deal with. But even with the reduction trick, it is difficult to surpass, 
say n = 7 for c = 2. Higher values of the coordination implies smaller values for the accessible height. 

Thus we conclude that exploring the behavior of T, for large n, dealing with exact results, is impossible. For this 
reason we have turned our interest towards developing simple approximate methods which can, however, provide 
interesting information on the asymptotic value of T. In this context appears the idea of obtaining bounds for T. It 
is found that by performing adequate averages of the decay widths appearing at each stage of breaking of a height n, 
the value of T obtained in the next height n -I- 1 is systematically lower (or higher) than what the exact result would 
be. In one Appendix we have given details of why bounds result. As the results obtained from the bounds reach 
values beyond n = 17 (c = 2), one is able to explore their asymptotic behavior by a careful exponential fitting which 
provides clear numerical evidence (although non rigorous) that for c = 2, T tends to a non-zero value when n tends 
to infinity. This conclusion is obtained for both the power-law and the exponential breaking rules. For the power 
law hazard rate, a proof was given in p9[ . Invoking conventional universality- class arguments one deduces that this 
non-zero limit holds for hierarchical structures of any coordination. 
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APPENDIX A: ON THE REPLICAS 

We have seen in Section III that the knowledge of the reduced configuration of a level n — 1 is not enough to obtain 
the primary diagram of the level n; we also have to know the replica of the n — 1 level. Expressed in the singular, 
this sentence is misleading. In fact it only holds for c = 2. One can easily check that for a general c , the number of 
replicas required, m, is 

1 <m< (c- 1) . (Al) 

In the case of the power-law breaking rule, any decay width of these replicas would be obtained by multiplying its 
normal value by the factor 

10 



(A2) 



while as seen in Section V, the exponential breaking rule demands the individualized calculations of each replica, with 
its corresponding extra loading. As an example beyond the usual c = 2, let us consider for the case c = 4 the process 
of breaking up to the collapse of the two minimum trees n — and n ^ I. Using a self-explanatory notation, we have: 



n = |T| 
n=l 

00 _ 

00 



0H 

as 



0J|_0 

2ir2 



0J|_0 

olH 



0J|_0 

olfo 



We see that the solution of the height n = 1, demands the information of the decay width of |jj but also that of 
of 2 and of 4 ; i.e. for c = 4 the iterative method requires the knowledge of three replicas, as foreseen in Eq. (All) 



APPENDIX B: ON WHY BOUNDS RESULT 



Following the arguments of section III and IV, one easily sees that the first S in which there must be a discrepancy 
between the exact result and the approximate results coming from the use of effective diagrams is the S^ oi n — 4. 
For c = 2, p = 2, we obtain 



(^3 (exact) 



128 

2893 

17 

385 " 



0.0442447 



0.0441558 






1 

IT 



1 



59 



1342 



8 + 202/5.143/5 
0.0439642 



= 0.0441074 



(Bl) 



This results from the fusion of the two configurations | 3.1 1 and 3.2 of Fig. 1^, which have a different F. To clarify 
why bounds result, let us analyze this point from a general perspective. In Fig. Wm is drawn the top of a reduced 
diagram of height, say, n, and at its right the corresponding effective diagram. We assume Fi ^ F2, 02 = 71 + 72, 
and i = AM, GM or HM, 



0-3, AM 



71 



ri 



72 1 



a2 a2 

03, GM = r^"" • F2" , 

03, AM 



(B2) 



1 



(il) J_ + (22.] J_ 

\a2j Ti ^ \a2j r2 



From the reduced diagram we obtain the top of the corresponding primary of the height n + I. This is shown in 
Fig. n2, and from the effective diagram one obtains the top of the primary shown in Fig. p^. Now let us compute the 
exact 63 coming from Fig. O to be compared with that (approximate) coming from Fig. |13|. In Fig. |l^, we have 

p(3.1) 

p(3.2) 

p(3.3) 



ai 


71 


ao + ai 


ao + a2 


ai 


72 


ao + ai 


a2 + a2 


Go 


ai 



ao 



ao + ai ao + ai a2 + a2 



Thus 



53(exact) = p(3.1] 



1 



ao + Ti 



p(3.2) 



ao 



p(3.3) 



1 



ai + 02 



(B3) 
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Analogously, in Fig. |3[ 



p(3.1) 
p(3.2) 



ai 02 






flo + ai ao + fli a2 + a2 



and 



53,,:=p(3.1) \ +p(3.2)- 



a3,i + flo 



fll + 02 



As the third term of Eq. (B3) coincides with the second of Eq. (B4) let us reorder Eq. (B3), giving 

(53 exact) -p 3.3 ■ = ■ ■ -— + -— 

ai + a2 (ao + ai)(ao + a2) Voo + ri ao + T 



and similarly in Eq. (B4), obtaining 



S3,^-pi2.2)- 



1 



ai • 02 



1 



oi + 02 (oo + ai)(oo + 02) Vo3,4 + oo 
To simplify the comparison, let us define two new functions, 



(B4) 



A3(exact) = I 53(exact) — p(3.3) 
A3(z)= (63. — p{2.2 



1 



(ao + l)(oo + 02) 7i 



72 



oi + a2y ai 

1 \ (ao + l)(ao + 02) 



Oo + Fi ao + F2 



Oi + 02 



02 



Ol 



Oo + 03,i 



(B5) 
(B6) 



Thus, the As represent the Ss after adding an equal term, and multiplied by an equal factor. It is interesting to observe 
the effect produced by the fusion of | 3.1 1 and 3.2 from an algebraic point of view: the sum of the two fractions of 



Eq. (B5) has converted into the fraction at the right of Eq. (B6). In the case of Fi — F2, A3(exact)= A3(i). In other 
words, if Fi = F2 that fusion is exact. 

To deal with the Arithmetic Mean, let us define 



/(^) 



02 



ao + X 



(B7) 



(oo, a2 > 0) which is concave; this implies that /(AiFi + A2F2) < Ai /(Fi) + A2 /(F2) where Ai = p-, A2 



— , Ai + A2 = 1 and therefore 



02 



Oo 



71 1 1 72 1 

ti2 Ti a2 r2 



< 



71 



72 



Oo + Fi Oo + F2 ' 



which means that A3(AAf) < A3 (exact) and therefore ^zi^AM) < ^3 (exact). 

To deal with the Geometric Mean let us in Eq. (B7) make the change of variable z = logx, then 



02 



02 + a; 
which is also a concave function in z. Hence we have 



5(^) 



02 



ao + e^ 



Ai- 



02 



02 + e^ 



+ A2 



02 



> 



02 



02 + e^2 - ao + e(Pi^i+P2Z2) ' 



71 



72 



> 



02 



02 



00 + Fi a2 + F2 - ao + (ei°g^i)^i (eiog^2)A2 ao + F/^F/ 



Thus A3(GA/) < A3 (exact) and S^iGM) < (53 (exact) 



Finally, for the higher bound we will make in Eq. (B7) the change of variable z = - : 



12 



/(x) = ■ — = h[z) 



02 + a; a^z + 1 

h{z) is a convex function in 2:, therefore 

/i(AiZi + A2Z2) > Ai ft-(zi) + A2 ^1(2:2) , 
and in terms of our ordinary variables this iniphes 

A3(if.M.) = / f pT^j > Pi /(ri) +P2 /(r2) = A3(cxact) 

and hence 63(11. M.) > (^3 (exact). 

Note that the argument presented is vaUd for the S3 of any n of coordination c — 2, and for both the power-law 
breaking rule and for the exponential breaking rule. 
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TABLE I. Main asymptotic results for the three standard modahties of FBM in the static and dynamic cases. 



ELS 



LLS 



HLS 



Static 



Dynamic 



CRITICAL POINT 

a, = e-i/" 

Daniels (1945) 

CRITICAL POINT 

T. = 1/p 

Coleman (1958) 



NO CRITICAL POINT 

ac oc 1/lnN 

Harlow (1985) 

NO CRITICAL POINT 

Kuo and Phoenix (1987) 



NO CRITICAL POINT 

(7c cc 1/lnlnN 

Newman and Gabrielov (1991) 

CRITICAL POINT 

This paper 



TABLE II. Sensitivity analysis of the exponential decay fitting to the lower bound results. 



20 



k 

"Tg" 

18 
17 
16 
15 
14 
13 
12 
11 
10 
9 
8 
7 
6 
5 



T, 



l.oc 



19 



T, 



Z.oo 



18 



T, 



l.oc 



0.32575 
0.32551 
0.32541 
0.32537 
0.32536 
0.32536 
0.32536 
0.32537 
0.32537 
0.32537 
0.32537 
0.32537 
0.32537 
0.32537 
0.32537 



18 
17 
16 
15 
14 
13 
12 
11 
10 
9 
8 
7 
6 
5 



0.32579 
0.32553 
0.32542 
0.32537 
0.32536 
0.32536 
0.32536 
0.32537 
0.32537 
0.32537 
0.32537 
0.32537 
0.32537 
0.32537 



17 
16 
15 
14 
13 
12 
11 
10 
9 
8 
7 
6 
5 



0.32585 
0.32555 
0.32542 
0.32538 
0.32536 
0.32536 
0.32536 
0.32536 
0.32537 
0.32537 
0.32537 
0.32537 
0.32537 
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FIG. 1. Comparison of Monte Carlo simulations; the broad distributions come from using the standard approach, and the 
thinner distributions come from using the probabilistic approach. The time to failure is plotted in dimensionless units. 

FIG. 2. Breaking process for the three smallest trees of coordination c = 2 {N — 1, 2, 4). The integers in parenthesis (r) 
represent the number of breakings occurred. The 5s stand for the time elapsed between successive individual breakings and the 
numbers under the legs indicate the load they bear. 

FIG. 3. Calculation of three partial decay widths in n = 3, from the information obtained in n = 2. 

FIG. 4. Symbolic representation of the gradual rupture of the tree of height n = 3 (c = 2, p = 2). Time flows downwards. 

FIG. 5. Result of "reducing" the diagram of Fig. U. 

FIG. 6. "Effective" diagram for n = 3. 

FIG. 7. Dimensionless lifetime, T, for a fractal tree of height n. The small circles are obtained from Monte Carlo simulations. 
Lines 4 and 1 are higher bounds based on Tmin and the HM respectively. Lines 2, 3 and 5 are lower bounds based on the GM, 
AM and Tmax respectively. 

FIG. 8. Visualization of the exponential fittings to the results obtained by using Geometric and Harmonic means. 

FIG. 9. Primary diagrams for the exponential breaking rule (n = 2). 

FIG. 10. Results for T from trees of height n, with the exponential breaking rule. The small circles correspond to Monte 
Carlo simulations. Lines 2 and 3 correspond to lower bounds based on GM and AM respectively. 

FIG. 11. Top of a failure diagram down to the fourth decay stage, (a) represents the reduced diagram and (b) the corre- 
sponding effective diagram. 

FIG. 12. Top of the primary diagram built from the reduced diagram drawn in Fig. hll (a). 

FIG. 13. Top of the primary diagram built from the effective diagram drawn in Fig. hll (b). 
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